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INTERNAL  ENERGY  MODE  RELAXATION  IN  HIGH 
SPEED  CONTINUUM  AND  RAREFIED  FLOWS 


By  Eswar  Josyula*  and  Dean  C,  Wadsworth* 


Introduction 

The  presence  of  shock  waves  in  high  speed  flow,  of  a  polyatomic  gas  presents  consid¬ 
erable  difficulties  for  accurate  numerical  simulation  of  the  flowfield.  The  shock  wave 
redistributes  the  high  kinetic  energy  of  the  oncoming  flow  into  various  internal  en¬ 
ergy  modes,  which  relax  relatively  slowly,  leading  to  significant  chemical  and  thermal 
nonequilibrium  in  the  stagnation  region.  In  the  gas  kinetic  description,  intermolecular 
collisions  change  the  translational,  rotational,  vibrational,  and  electronic  energies  of 
the  collision  partners.  The  probabilities  or  effective  cross  sections  of  these  elementary 
processes  differ  significantly,  giving  rise  to  widely  separate  relaxation  times  for  the 
internal  modes.  Thus  it  becomes  important  to  account  for  the  rates  of  relaxation  pro¬ 
cesses  to  predict  the  nonequilibrium  behavior  of  these  kinds  of  flows.  The  continuum 
description  is  well  suited  at  lower  altitudes  of  the  flight  regime  for  the  prediction  of 
aerodynamic  loads  and  heating  rates  on  the  thermal  protection  systems.  However, 
at  high  altitudes  and  associated  low  densities  the  larger  mean  free  path  Invalidates 
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the  continuum  assumption  and  the  rarefied  solution  approaches  are  necessary  [1], 
Among  the  solution  approaches  in  hypersonic  rarefied  flows,  the  Direct  Simulation 
Monte  Carlo  (DSMC)  method  is  widely  used  [2]. 

In  the  continuum  approach  the  Navier-Stokes  equations  consist  of  reaction  prob¬ 
abilities- for  quantifying  the  thermal  and  chemical  nonequilibrium  effects.  These  are 
typically  temperature-dependent,  though  the  concept  of  temperature  in  a  nonequilib¬ 
rium  flow  is  ill-defined.  In  the  DSMC  method,  however,  the  internal  energy  relaxation 
occurs  as  a  consequence  of  intermolecular  collisions.  This  requires  knowledge  of  the 
level-to-level  cross-sections  for  internal  energy  modes  and  the  energy-dependent  cross- 
sections  of  chemical  reactions.  Recent  studies  [3, 4]  were  aimed  at  finding  a  consistent 
way  of  interpreting  nonequilibrium  in  both  the  approaches.  The  work  of  Ref.  [4]  de¬ 
veloped  an  exact  relationship  that  connects  the  vibrational  relaxation  number  used  in 
DSMC  method  to  the  continuum  method,  A  meaningful  comparison  of  the  reaction 
probabilities  and  cross-sections  is  critical  to  the  development  of  numerical  predictive 
tools  that  span  the  continuum  and  rarefied  regimes.  It  is  one  of  the  objectives  of  the 
proposed  study. 

The  internal  energy  mode  relaxation  in  the  shock  wave  structure  was  studied  using 
continuum  and  rarefied  approaches  by  various  investigators.  See  for  example,  Refs.  [5, 
6].  The  advantage  of  the  shock  structure  problems  is  that  it  does  not  involve  the 
nonequilibrium  effects  of  gas-surface  interactions,  an  area  of  research  not  considered 
in  the  proposed  paper.  Numerical  experiments  performed  by  Pham-Van-Diep,  et 
al.  [7]  in  testing  continuum  descriptions  showed  that  the  continuum  solutions  of  the 
shock  structure  in  a  monatomic  gas  were  able  to  match  the  rarefied  solution  at  a 
Mach  number  of  1.2.  At  higher  Mach  numbers  the  solutions  obtained  with  continuum 
approach  do  not  match  those  by  the  DSMC  method.  The  proposed  study  will  address 
the  inconsistencies  in  relaxation  mechanisms  implemented  in  Continuum  and  DSMC 
codes. 

The  rotational  energy  relaxation  in  shock  structures  of  Nitrogen  was  studied  using 
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Burnett  equations  by-  Lumpkin,  et  al.  [8],  The  model  was  validated  with  experiment 
up  to  Mach  6.  In  the  vibrational  energy  relaxation  implemented  in  the  continuum  ap¬ 
proach  of  Josyula.  [9]  the  master  equations  were  coupled  to  the  fluid  dynamic  equations 
to  model  the  nonequilibrium  flow  physios  in  Hypersonic  flow.  These  two  modeling 
approaches  will  be  used  in  the  present  study  to  compare  with  results  obtained  by  the 
DSMC  code. 


Analysis 


This  section  gives  the  governing  Navier-Stokes  equations  used  for  coupling  to  the 
Lumped  Landau-Teller  Vibrational  Relaxation  Model  (LLTR)  and  Discrete  State 
Kinetic  Relaxation  Model  (DSKR). 

The  global  conservation  equations  in  mass-averaged  velocity  form  are  shown  be¬ 
low: 


+  V  ■  [pt,(u  +  Un2  +  uvn2)]  =  d>v  v=0,l,...  (1) 

Q 

-Qt{Ps)  +  V  •  [ps(u  +  US)]  =  U!s  (2) 

d 

-(pu) +  V- (pirn +  r)  =  0  (3) 

d 

g^{pevib)  ^  Us)  H-  Qvih]  —  P&vib  &vib^  Qt—v  (4) 


d 

—  (pe)  +  V  •  [p(e  +  p/p) U  -  (Y,  ()mb  +  qtrans)  +  £(M.U»)  ~  U  ■  f]  =  0  (5) 

The  conservation  Equation  (1)  used  in  the  DSKR  code  is  written  for  mass  density  in 
quantum  level  v  for  diatomic  nitrogen.  The  term  uN2  denotes  the  diffusion  velocity  of 
component  N2  of  the  gas  mixture  and  uvn2  is  the  diffusion  velocity  of  level  v  relative 
to  N2  diffusion  velocity.  The  source  term  uv  derived  from  the  vibrational  master 
equations  is  made  up  of  the  relevant  energy  exchange  processes  consisting  of  the  V-T 
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and  V-V  reaction  mechanisms.  Of  the  three  species  (02,  Ar2,  and  0)  considered  for 
the  air  mixture  in  the  DSKR  code,  only  the  species  N2  was  treated  as  an  anharmonic 
oscillator  in  the  DSKR  model  with  the  following  energy  exchange  mechanisms. 

—  (Pn2- N2Y  T+(pN2-N2y  1  +(PW2-02)'  r+(pW2-o)U  T +  (PJV2->Ar)*SS  +  (pN->N2YeCOmb 

(6) 

The  density  of  molecular  nitrogen  is  the  sum  of  population  densities  in  the  various 
vibrational  levels. 

Pn2  =  Y,  Pv  (7) 

»=o,i.. 

The  mass  conservation  of  species  treated  in  the  LLTR  model  is  represented  by  Equa¬ 
tion  (2).  The  production  of  small  amounts  of  atoms  due  to  dissociation  of  molecules 
is  included  in  the  source  term,  ws.  The  term  us  denotes  the  diffusion  velocity  of 
component  s  of  the  gas  mixture.  The  mixture  density,  p  is  the  sum  of  the  partial 
species  densities, 

P  =  Pn2  +  Po2  +  Po  +  Pno  +  Pn  (8) 

For  simulations  in  air,  the  maximum  temperature  of  the  shock  was  below  7,000  K, 
thereby  necessitating  the  oxygen  dissociation  reaction  in  the  DSKR  code.  Due  to 
small  amount  of  nitrogen  dissociation  and  the  difficulty  of  incorporating  collision 
probabilities  of  additional  species,  only  oxygen  dissociation  was  considered  in  the 
DSKR  model.  However,  the  LLTR  code  considered  N  and  NO  also. 

Equation  (3)  gives  the  conservation  of  total  momentum.  Equation  (4)  is  the  con¬ 
servation  equation  for  vibrational  energy  where  Qt-v  denotes  the  energy  exchange 
between  the  vibrational  and  translational  modes.  For  diatomic  nitrogen  in  the  DSKR 
code,  a  separate  vibrational  conservation  equation  was  not  necessary  as  the  vibra¬ 
tional  energy  was  calculated  at  each  quantum  level,  discussed  later.  The  diatomic 
oxygen  molecules  in  the  DSKR  code  and  all  the  species  in  the  LLTR  model  were 
assumed  as  harmonic  oscillators.  Equation  (4)  was  solved  with  a  source  term  for  the 
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V-T  coupling,  modeled  according  to  the  Landau-Teller  [10,  11]  form: 


Vyib 


evibs  eviba 


where  evibs  is  the  vibrational  energy  of  the  molecular  species  s  and  e*ibs  is  the  vi¬ 
brational  energy  in  thermal  equilibrium  at  the  local  translational  temperature,  the 
relaxation  time  given  by 


Es  Xs/tlt 


(10) 


where  Xs  denotes  the  species  mole  fraction  and  tlt  is  the  Landau-Teller  inter-species 
relaxation  times.  tLt  was  computed  using  the  expression  developed  by  Millikan  and 
White  [12].  The  vibrational  temperature  of  molecular  species  s  was  determined  by 
inverting  the  expression  for  the  vibrational  energy  contained  in  a  harmonic  oscillator 
at  the  temperature,  Tv, 


6vib{s)  ~  ee„(S)/T,.  _  i 


(11) 


where  R  is  the  species  gas  constant  per  unit  mass.  V-V  exchanges  were  considered 
only  between  the  nitrogen  molecules  in  the  DSKR  model  used  to  compare  to  the 
LLTR  model;  they  were  neglected  for  the  oxygen  molecules  in  the  DSKR  model  and 
all  the  species  in  the  LLTR  model.  The  DSKR  model  incorporating  the  new  vibration- 
dissociation  coupling  neglected  V-V  exchanges  also.  Equation  (4)  also  includes  terms 
for  the  conduction  and  diffusion  of  vibrational  energy.  The  conservation  of  total 
energy  is  given  by  Equation  (5)  with  heat  conduction  and  species  diffusion  terms. 


The  kinetics  of  the  particle  exchanges  among  the  quantum  states  of  V2  were 
simulated  by  the  vibrational  master  equations.  The  population  distributions  were 
calculated  by:  [9] 


U)v 


=  [kvAv1  ->  v)pV'P  -  kvr(v  -+  v')pvp] 

+  ^2  [kvv(v',  w'  — >  v.  w)pv>pw>  -  kVv{v,  w  v ',  w')pvpw]  > 

w  j 


(12) 


where  only  single  quantum  transitions  have  been  considered.  The  equations  governing 
the  V-T  reactions  responsible  for  the  variation  of  the  particles  distributed  in  the  vth 


vibrational  level  of  diatomic  nitrogen  are: 


N2(v)  +  M  ^  N2(v  -  1)  +  M  (13) 

where  M  represents  02,  O,  N2.  The  equations  governing  the  V-V  processes  in  N2 
giving  the  reactions  responsible  for  the  variation  of  the  particles  distributed  in  the 
vth  vibrational  level  are: 


N2(v)  +  N2(w )  N2(v  —  1)  +  N2(w  +  1)  (14) 

For  the  kinetics  of  diatomic  nitrogen,  the  present  study  used:  (a)  V-T  forward  rate 
coefficients  calculated  according  to  expressions  proposed  by  Capitelli,  et  al,  [13]  and 
Billing  and  Fisher  [14]  and  (b)  V-V  forward  rates  by  Doroshenko,  et  al.  [15],  The 
V-T  forward  rate  coefficient  for  N2-0  collisions  was  from  the  work  of  Capitelli  [16] 
which  was  based  on  Refs.  [17]  and  [18].  Reverse  rate  coefficients  were  derived  from 
detailed  balance. 


The  vibrational  energy  of  the  N2  molecule  is  given  in  terms  of  the  quantum  level 
energies  by 

Pi 


evibN2  =  ~ei 

i= 1,2,..  PN 2 


(15) 


where  the  index  i  is  used  to  denote  the  quantum  level.  In  this  equation,  is  the 

pn2 

fractional  population  of  the  i vibrational  level  and  e^  the  quantum  level  energy  given 
by  the  third-order  approximating  formula: 


he  2)  i— 1,2,.. .1+1 


(16) 


The  above  equation  represents  anharmonic-oscillator  behavior  of  the  N2  molecule, 
where  h  is  the  Planck’s  constant  and  c  is  the  speed  of  light.  The  spectroscopic 
constants  are  given  by,  [19]  we= 2358.57  crn~\  wea;e=14,324  cm~ \  and  weye=-0. 00226 

cm-1.  When  45.  the  value  of  energy  exceeds  the  N2  dissociation  energy,  9.62 
eV.  [20] 
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Numerical  Procedure 


In  the  continuum  approach,  the  Roe  approximate  Riemann  solver  is  implemented 
in  finite  volume  formulation  by  computing  the  cell  interface  flux  as  a  summation  of 
wave  speeds  as  described  by  Cinnella,  et  al.  [21].  The  second  order  spatial  accuracy  is 
obtained  by  employing  the  MUSCL  approach  in  conjunction  with  the  minmod  limiter 
to  reduce  the  solution  to  first  order  accuracy  in  the  vicinity  of  strong  shock  waves, 
as  described  in  the  work  of  Josyula,  et  al  [22],  The  entropy  correction  for  the  Roe 
scheme  is  implemented  as  discussed  in  the  Reference  [22],  The  viscous  terms  are 
evaluated  using  central  differencing.  An  explicit  predictor-corrector  method  is  used 
to  advance  the  solution  in  time.  This  approach  was  discussed  for  the  flux- splitting 
option  by  MacCormack  in  Reference  [23], 

The  rarefied  approach  uses  the  Direct  Simulation  Monte  Carlo  (DSMC)  code 
developed  by  Wadsworth,  et  al.  [24],  The  code  was  used  to  study  nonequilibrium 
flowfield  applications  and  includes  energy  exchange  mechanism  for  vibrational  and 
rotational  energy  transfers. 


Preliminary  Results  and  Discussion 

The  Navier-Stokes  continuum  code,  described  above,  using  the  perfect  gas  assumption 
was  used  to  compute  the  internal  structure  of  a  shock  wave.  Comparisons  were 
made  for  low  Mach  Number  shock  waves  with  (a)  theoretically  obtained  Navier- 
Stokes  solution,  (b)  available  experimental  data  and  (c)  the  results  obtained  by  a 
DSMC  code.  These  perfect  gas  computations  and  comparisons  were  made  to  establish 
a  baseline  for  which  the  internal  energy  modes  of  rotational  and  vibrational  mode 
relaxation  will  be  computed  in  the  proposed  paper. 

Figure  1  shows  the  comparison  of  flow  parameters  inside  a  Mach  2  shock  wave 
using  a  perfect  gas  Navier-Stokes  computational  code  and  the  exact  Navier-Stokes 
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solution.  For  the  perfect  gas  assumption,  the  agreement  of  both  sets  of  results  is 
excellent. 

Experimental  measurements  reported  by  Sherman  [25]  inside  the  shock  wave  struc¬ 
ture  were  used  for  comparison  in  the  next  figure.  Temperature  comparisons  using  the 
present  Navier-Stokes  computational  code  using  the  perfect  gas  assumption  and  the 
experimental  data  are  shown  in  Fig.  2.  The  flow  Mach  number  is  1.9  in  air  and  the 
computed  results  are  shown  for  three  different  bulk  viscosity  coefficients,  A  =  -2/3/i. 
A  =  +4/3 n,  and  A  =  2/i.  The  internal  energy  relaxation  related  to  the  bulk  mod¬ 
ulus  does  have  a  significant  effect  on  the  result.  For  the  commonly  used  value  of 
A  =  -2/3 n,  the  temperature  is  underpredicted  for  M  <  1  and  overpredicted  for 
M  >  1.  Increasing  the  value  of  A,  however,  compensates  for  this  under-  and  over¬ 
prediction.  This  signifies  the  importance  of  rotational  and  vibrational  relaxation  in 
the  shock  structure. 

Figures  3  to  7  show  the  variation  of  pressure,  density,  temperature  and  velocity 
Argon  shock  structure  using  the  perfect  gas,  Navier-Stokes  code  and  DSMC  method. 
The  continuum  calculations  for  Mach  1.2  (Fig.  3)  match  those  obtained  using  the 
DSMC  method.  However,  at  a  higher  Mach  number  of  2  shown  in  Fig.  4,  one  can  see 
the  thicker  shock  profile  for  the  pressure  variation  across  the  shock  wave  obtained  by 
the  DSMC  code  compared  to  the  Navier-Stokes  code.  The  mass  density  obtained  by 
the  DSMC,  Fig.  5,  shows  the  thicker  shock  thickness,  the  thickness  slightly  higher  for 
M  >  1.  The  temperature  and  velocity  profiles  of  the  continuum  solutions,  Fig.  6  and 
Fig.  7  show  greater  thickness  in  the  region  for  M  <  1. 


8 


References 


[1]  E.P.  Muntz.  Rarefied  Gas  Dynamics.  Annual  Review  of  Fluid  Mechanics , 
21:387-417,  1989. 

[2]  G.A.  Bird.  Molecular  Gas  Dynamics  and  the  Direct  Simulation  of  Gas  Flows. 
Oxford  University  Press,  inc.,  New  York,  1994. 

[3]  M.S.  Ivanov  and  S.F.  Girnelsliein.  Computational  Hypersonic  Rarefied  Plows. 
Annual  Review  of  Fluid  Mechanics.  30:469-505,  1998, 

[4]  N.E.  Gimelshein,  S.F.  Gimelshein,  and  D.  A.  Levin.  Vibrational  Relaxation  Rates 
in  the  Direct  Simulation  Monte  Carlo  Method.  Physics  of  Fluids,  14(12):4452- 
4455,  December  2002. 

[5]  D.  Gilbarg  and  D.  Paolucci.  The  Structure  of  Shock  Waves  in  the  Continuum 
Theory  of  Fluids.  Journal  of  Rational  Mechanics  and  Analysis,  2(4):617-642, 
1953. 

[6]  G.A.  Bird.  The  Velocity  Distribution  Function  Within  a  Shock  Wave.  Journal 
of  Fluid  Mechanics,  30(3):479-487,  1967. 

[7]  G.C.  Pham-Van-Diep,  D.A.  Erwin,  and  E.P.  Muntz.  Testing  Continuum  De¬ 
scriptions  of  Low  Mach  Number  Shock  Structures.  Journal  of  Fluid  Mechanics, 
232:403-413,  1991. 

[8]  F.E.  Lumpkin  III,  D.R.  Chapman,  and  C.  Park.  A  new  rotational  relaxation 
model  for  use  in  hypersonic  computational  fluid  mechanics.  AIAA  Paper  89- 
1737,  1989. 

[9]  E.  Josyula.  Computational  Study  of  Vibrationally  Relaxing  Gas  Past  Blunt  Body 
in  Hypersonic  Flows.  Journal  of  Thermophysics  and  Heat  Transfer,  14(1):18— 26, 
January  2000. 


9 


[10]  L.  Landau  and  E.  Teller.  Zur  Theorie  der  Schalldispersion.  Physikalische 
Zeitschrift  der  Sowjetunion,  10(l):34-43,  1936. 

[11]  W.G.  Vincenti  and  C.H.  Kruger  Jr.  Introdxiction  to  Physical  Gas  Dynamics. 
John  Wiley  h  Sons,  1967. 

[12]  R.C.  Millikan  and  D.R.  White.  Systematics  of  Vibrational  Relaxation.  Journal 
of  Chemical  Physics ,  39(12):3209-3213,  1963. 

[13]  M.  Capitelli,  C.  Gorse,  and  G.D.  Billing.  V-V  Pumping  Up  in  Nonequilibrium 
Nitrogen:  Effects  on  the  Dissociation  Rate.  Chemical  Physics ,  52(3):299-304, 
1980. 

[14]  G.D.  Billing  and  E.R.  Fisher.  VV  and  VT  Rate  Coefficients  in  Diatomic  Nitrogen 
'  by  a  Quantum  Classical  Model.  Chemical  Physics ,  43(3) : 395—40 1 ,  1979. 

[15]  V.M.  Doroshenko,  N.N.  Kudryavtsev,  S.S.  Novikov,  and  V.V.  Smetanin.  Effect 
of  the  Formation  of  Vibrationally  Excited  Nitrogen  Molecules  in  Atomic  Recom¬ 
bination  in  a  Boundary  Layer  on  the  Heat  Transfer.  High  Temperature  (USSR), 
28(l):82-89,  1990. 

[16]  M.  Capitelli,  I.  Armenise,  and  C.  Gorse.  State-to-State  Approach  in  the  Kinetics 
of  Air  Components  Under  Re-Entry  Conditions.  Journal  of  Thermophysics  and 
Heat  Transfer,  ll(4):570-578,  October  1997. 

[17]  G.D.  Billing.  VV  and  VT  Rates  in  N2-02  Collisions.  Chemical  Physics, 
179(3):463-467,  1994. 

[18]  A.K.  Mnatsakanyan  and  G.V,  Naidis.  The  Vibrational-Energy  Balance  in  a 
Discharge  in  Air.  High  Temperature  (USSR),  4(3) :506— 513,  1985. 

[19]  K.P.  Huber  and  G.  Herzberg.  Constants  of  Diatomic  Molecules.  Van  Nostrand 
Reinhold  Company,  New  York,  1979. 


10 


V.  |l 


[20]  I.  Armenise,  M.  Capitelli,  R.  Celiberto,  G.  Colonna,  and  C.  Gorse.  The  Effect  of 
N  +  N2  Collisions  on  the  Non-equilibrium  Vibrational  Distributions  of  Nitrogen 
under  Reentry  Conditions.  Chemical  Physics  Letters ,  227(1):  157— 163,  1994. 

[21]  P.  Cinnella  and  B.  Grossman.  Upwind  Techniques  for  Plows  with  Multiple  Trans¬ 
lational  Temperatures  .  AIA A  Paper  90-1660,  June  1990. 

[22]  E.  Josyula.  and  J.S.  Shang.  Computation  of  Hypersonic  Flowfields  in  Thermal 
and  Chemical  Nonequilibrium.  AIAA  Paper  92-2874,  July  1992. 

[23]  R.  MacCormack.  Current  Status  of  Numerical  Solutions  of  the  Navier-Stokes 
Equations.  AIAA  Paper  85-0032,  1985. 

[24]  D.C.  Wadsworth,  I.J.  Wysong,  and  D.H.  Campbell.  Direct  Simulation  Monte 
Carlo  Modeling  for  Non-Equilibrium  Flowfield  Applications.  Technical  Report 
PL-TR-97-3024,  Air  Force  Research  Laboratory,  Edwards  Air  Force  Base,  CA, 
May  1999. 

[25]  F.S.  Sherman.  A  low  density  wind-tunnel  study  of  shock- wave  structure  and 
relaxation  phenomena  in  gases.  Technical  Report  Technical  Note  3298,  National 
Advisory  Committee  for  Aeronautics,  Washington,  July  1955. 


11 


Flow  variable 


Figure  1:  Comparison  of  flow  parameters  across  shock  wave  in  Mach  2  air  flow  between 
Navier-Stokes  CFD  calculation  and  exact  solution 
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Sherman’s  Air  experiments 
translational  temperature 


X=  4.28  x  10'4  m 


Figure  2:  Comparison  of  flow  parameters  across  shock  wave  in  Mach  1.9  air  flow 
between  Navier-Stokes,  perfect  gas  CFD  solution  and  Experiment  showing  effect  of 
varying  the  bulk  viscosity  in  the  Navier-Stokes  equations 
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Flow  variable 


Figure  3:  Variation  of  flow  parameters  across  shock  wave  in  Mach  1.2  Argon 
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Argon  gas  static  pressure 


X  =  4.0261 49  xl  O'6  m 


Figure  4:  Comparison  of  static  pressure  across  shock  wave  in  Mach  2  argon 
between  Navier-Stokes  CFD  solution  and  DSMC 


Figure  5:  Comparison  of  mass  density  across  shock  wave  in  Mach  2  argon  flow  between 
Navier-Stokes  CFD  solution  and  DSMC 


Figuie  6.  Comparison  of  temperature  across  shock  wave  in  Mach  2  argon  flow  between 
Navier-Stokes  CFD  solution  and  DSMC 


